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ABSTRACT 

This  revised  version  of  CAC  Document  #21  supercedes  the  document  dated 
November  8,  1971. 

Several  methods  have  been  proposed  to  enable  the  computation  of  eigen- 
values and  eigenvectors  of  large,  real  symmetric  or  complex  Hermitian 
matrices  on  ILLIAC  IV. 

One  of  the  most  effective  methods  in  the  utilization  of  parallel  compu- 
tations has  proven  to  be  a  modified  Jacobi  algorithm.   This  document 
presents  yet  another  modification  which  exploits  the  parallelism  of 
ILLIAC  IV  to  a  greater  extent  than  has  been  previously  done. 

Flow  charts  and  the  assembly  language  (ASK)  routine  JACOBI  are  included 
in  the  report . 
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1.0   Introduction 

The  computation  of    eigenvalues   and   eigenvectors   of   large,  real   symmetric 
matrices     is    of  great   practical  value   in  many  fields.      One  of   the  most 
effective  methods   that  has   been   examined    to    solve   this   problem   on   the 
ILLIAC   IV   is  a  parallel  version  of   Jacobi's   algorithm. 

The  main  difference  between  this   code  and   the  one  previously  written   [1] 
is   that  within  each   sweep,   defined  by   (2m-l)    transformations    (where  m  = 
[(n+l)/2],    i.e.    the  greatest    integer    less    than  or   equal    to    (n+l)/2,    in 
which  n   is   the  order  of   the  matrix),    each  orthogonal   transformation 
annihilates  different    [-=-]    off -diagonal   elements.    This  proved   to  be  a 
substantial    improvement    that   led    to   a  greater   speed   of    convergence. 
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2.0       Jacobi' s  Method 


In  the  classical  method   of   Jacobi,    a  real   symmetric  matrix   is 


reduced    to    the  diagonal   form  by   a   sequence  of   plane  rotations 

A  ,       =  R,    A,    -    R  (k  =  1,2...),   where  A  =  A  is   the  original  matrix, 

and   the  rotation  matrix  R,    eliminates    the  off-diagonal   element  a 

(k)  (k)  Pq 

(hence  a  )through  an  angle    a  [5].      See  Appendix  A  for   the   appropriate 

pq  (k)  pq  (k) 

value  of  a  to   annihilate   the  element   a 

pq  pq 


3. 0     Modifications   to   the  Classical   Jacobi  Method 

3. 1  Decomposition   Into   Block  Diagonal   Submatrices 

Rather   than  searching   for   the  largest  off-diagonal   element  of  A 
in  the    (p,q)  position  and   eliminating   a        and  a      ,   A.    Sameh  and   L.    Han. [4] 

pq  qp 

proposed  a  modified  Jacobi  algorithm  where  all  off-diagonal  elements 
of  each  2x2  submatrix  along  the  diagonal  are  eliminated  through  an 
orthogonal   transformation. 

In  order   to  bring    to   the  diagonal   new  submatrices  with  non-zero   off- 
diagonal    elements,    the  method   necessitates   a   large  degree   of   row  and 
column  shuffling   which   tends    to   be   expensive  on  a   parallel   computer   such 
as    ILLIAC    rv. 

3. 2  Optimal    Construction  of   Orthogonal   Transformations 

A.  SamehP]  showed  that  a  judicious  choice  of  the  pairs  (p,q)  can 
produce  a  modified  Jacobi  algorithm  that  attains  maximum  efficiency  of 
parallel   computation. 

For   example,    for   a  matrix  A  of   order  4,    if    the  orthogonal   transfor- 
mation R   is   chosen  as, 


(3.1) 


R   = 


0 


'1 


0 


-s. 


-s, 


where  c^  =  cos  a^(i  =  1,2),  then  RAR1-  would  have  zero  elements  in  positions 
(1,3)  and  (2,4)  provided  that  the  angles  a]_  and  012  are  properly  chosen. 
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Define  m=    [(n+l)/2]   where  n   is    the   order   of    the  matrix  and    [    ]    is    the 

greatest    integer   function.      Let   each    of   the    (2m-l)    orthogonal   transformations 

2 
be  denoted   by   a   sweep.      Noting    that    the  maximum   number   of    the    (n   -n)/2 

off-diagonal    elements  which   can  be   eliminated   by   an  orthogonal    transfor- 
mation of    the   type    (3.1)    is    [n/2],    an  optimal   algorithm   requires    that,    [3]: 

1)  Each   orthogonal    transformation  R,     should    eliminate    [n/2]    off- 
diagonal    elements. 

2)  Each  sweep  should  eliminate  each  off-diagonal  element  once,  i.e. 
each  of  the  2m-l  orthogonal  transformations  in  a  sweep  should  annihilate 
different    [n/2]    off-diagonal    elements. 


3.2.1      Elimination   Scheme 

In    [3]    several    schemes  were  proposed    to   satisfy   the  requirements 
discussed    in   the   previous    section.      The   scheme    implemented    in   the  JACOBI 

algorithm  is   described  below: 

For   a  given  sweep,    each   of    the    (2m-l)    orthogonal  matrices   R, 
consist   of    the   elements, 

R(k)      =  R(k)=     _cos   a(k)      R(k)  = 

pp  qq  pq       pq         qp  pq 

(k) 

p    :•    q 


where   p   and   q    are   sequences   defined   by 

(a)  for    k   =   1,2,. . . ,m   -   1, 

q   =  m   -   k  +  1 ,    m-k+2,    ...  ,n-k. 

,,                  p=(2m-2k  +  l)-q,  m  -  k  +  1<_  q  <_  2m  -   2k, 

=    (4m    -    2k)    -    q,  2m    -    2k   <    q    ^  2m   -   k    -   1 , 

=    n,  2m    -   k   -    1    <    q , 

(b)  for    k   =  m,    m   +   1,...,    2m   -    1, 

q=4m-n-k,    4m-n-k+l,     ...     ,3m-k-l 

p  =  ii,  q  <  2m  -  k  +  1, 

(4m    -    2k)    -    q,  2m    -   k   +    I     :    q    ••    4m   -   2k   -   1 , 

=    (6m    -    2k    -    I)    -   q,  4m    -    2k    -    1    <    q . 


qp 

.         (k) 
=  sin  a 

pq 

.        (k) 

=  -sin  a 

pq 

-/.- 


(k) 
The  remaining  elements  of  R,  are  zero  except  for  n  odd,  then  R_   ,  _  .  =  1, 
°  k  2m-k,  2m-k 

For  a  given  k,  the  angles  a   are  determined  for  all  (p.q)  such  that  a 

pq  vtMH  pq 

(k) 
eliminates   the   element  a        ;    see  Appendix  A. 

pq  w 


For   example,    in  a  given  sweep,   denoting   each   element    in  the   transformation 
by   the   integer   k,    the  patterns  of   the  eliminated   elements   for  matrices   of 
orders   8   and   7  are  shown  below. 


3      6      2      5     14  17 

*     2      5     1     4      7  i  6 

*     1     4      7      3  ]  5 

*      7      3      6   14 

*      6      2   ,  3 

*      5   12 

*    i1 
!* 

7  X    llSy/ 

8  x  8i^ 


(k) 


(k) 


Note  that  since  av*    as  well  as  av"    is  eliminated,  if  one  completes  the 

qp  pq 

lower  diagonal  portion  of  the  matrix  above,  it  is  evident  that  any  given  k 
appears  in  each  row  and  column  once  and  only  once. 

For  further  examples  of  particular  orthogonal  transformations  constructed 
by  this  elimination  scheme  see  Appendix  B  . 
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4.0  JACOBI  -  An  ILLIAC  IV  Routine 

4. 1  Introduction 

JACOBI  is  an  ILLIAC  IV  routine  written  in  the  assembly  language  ASK, 
which  implements  the  modified  Jacob i  algorithm  discussed  in  Section  3.2. 

The  program  accepts  as  input  a  matrix  A,  and  n,  the  order  of  the 
matrix,  and  returns  as  output  the  matrix  A  reduced  to  diagonal  form,  the 
matrix  of  eigenvectors  corresponding  to  the  computed  eigenvalues  and  the 
number  of  sweeps  required  to  achieve  convergence. 

A  flow  chart,  a  description  of  JACOBI  and  auxiliary  routines,  and  a 
short  discussion  of  JACOBI  for  the  potential  user  are  now  presented. 

4. 2  Procedure 

Each  sweep  requires  2m-l  orthogonal  transformations  as  described  in 
Section  3.2.   For  each  transformation  the  following  sequences  of  events 

occur: 

(k-1) 
(i)      The  pairs    (p,q)    corresponding    to   the  element  a  to  be 

eliminated  are  determined. 

(ii)      The  orthogonal    transformation  matrix  R,     is    constructed    in 

order    to   eliminate   the  proper   elements  of   A. 

(iii)      A  is   pre-  and  post-multiplied  by  the   transformation  matrix 

t  (k) 

R,     to   yield   A,     =  R.       A,     -.    R.    with    elements   a  =0. 

k  k  k        k-1      k  pq 

(iv)      E,  _1  ,    the   eigenvector  matrix,    is   pre-multiplied  by  R,     to  yield 

E,  ,   where  E      is    the   identity  matrix.    (Note   that   the  rows  and   columns 
k  o 

of    E  correspond   to   the  left  and   right   eigenvectors   respectively) . 
After    2m-l   such    transformations   have  been  applied,    the   following 
convergence  criterion  are  subjected   to  A   : 

(a)  if    the  sum  of    the  squares   of    the  off-diagonal   elements    is 
zero,    convergence   is   attained. 

(b)  if    the  ratio   of    the  sum  of    the   squares  of    the  off-diagonal 

elements    to   the  sum  of    the  squares   of    the  diagonal   elements 

—  ~\(\ 
at   step     k    is   sufficiently  small    (10        )    in   comparison   to  an 

equivalent   ratio   at   step   1,    the  method  has   converged.      (to 

insure  numerical    stability   this   critera   is   not  applied 

for    the   first    three  sweeps). 
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If    the   convergence  tests   fail,    a  new  sweep   is   initiated.      When  the 
method   does   converge,    the  bounds   on   the   eigenvalues    are   computed 
using  Gerschgorin  discs, and   this    information  is   output  for   the  user 

4. 3     Main  Flowchart 

Notation:      A^   -  Matrix   to  be  diagonalized   at   step   k 

E,    -  Corresponding   eigenvector  matrix 

R,     -   Orthogonal    transformation  matrix 

B      -   Temporary  matrix 

Ratlo!     J.  (k)    2     /l       ak>2     where  a(k)    e  A, 


Kconv:      I  (1)    2 


a 


£  -  Transformation  count 
swp  -  Sweep  count 

k  -  swp  x(2m-l)  +  £    where  m  =  I   -  J 

n  =  order  of  matrix 

A  detailed  discussion  of  the  components  of  JACOB I,  depicted  in  the 
flowchart  presented  in  this  section  may  be  found  in  Section  A. 4. 
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fAo input 

matrix 
n  -   order 
of  matrix) 


(^   JACOB I  ) 


Initialize 
Set  E  =1 


m 


swp  =  0 
kconv  =  O.C 


Set   up  matrix  defining 
pattern  of  pairs    (p,q) 
to  be   eliminated  for 
each   sweep. 


£  =   1,2,. . .2m-l 
note   k  =   swp   x    (2m-l)    +   I 


swp=swp+l 


Retrieve    (p,q)   patterr 
in  fcth  TOV-  Qf  matrix 
computed  above. 


1 


ANGLE 

^Compu  t  e   s  ines  &N 

:osines    of   R   ty 


B    is    skewed    in    this 
routine    to    prepare 
for    transpose   below 


MULTPL 


—  —  —  _  /      Compute 


Ote  A 


k-l 


k-l 


Align   rows   of 

,  t 

B      to    produce 

B      Vi  Rk 


see    (1)    p.    10 


see    (2)    p.    10 


see    (3)    p.    11 


see    (3)    p.    12 


0 
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Determine 
new  An 


Update  Eigenvector 
matrix. 


End   of    sweep    test 


Insure  numerical 
stability 


Convergence  Test 


see  (4)  p-  12 


see  (5)  P.  13 


X 


Find  Convergence 
Ratio 


see    (6)     P-  13 


GERSH 


Compute  Bounds 

on 
.Eigenvalues 


see    (7)      p.  13 


f    RETURN  j 
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4.4  Description  of  Main  and  Auxiliary  Routines 

After  saving  the  necessary  registers,  return  addresses,  and  addresses 

of  calling  parameters,  setting  up  constants  and  counters,  JACOBI  constructs 

E  =1,  the  identity  matrix, 
o 

(1)   Determination  of  pairs  (p,q) 

In  order  to  implement  this  phase  of  the  algorithm  on  ILLIAC 
IV,  it  is  desirable  to  maintain  compatibility  with  PE  numbering 

(i.e.  0<p,q<n-l)  and  also  alter  the  definition  for  p,q.   The 
following  definition  will  produce  identical  pairs  except  that  now 


a  .  . 
section. 


in  Section  3.2.1  corresponds  to  a..  ._.  as  defined  in  this 


Let  q  =  PE  number  q  =  0,1,... n-1 

i+ 
2 


m  =   [—5—]  n  =  dimension  of  matrix 


For  k  =  l,2,...2m-l    let  k  =  4m-2        k  =  l,2,...m-l 

o 

k  =  6m-3        k  =  m,m+l, . . .2m-l 
o 

Then  p  =  (k  -  2k  -  q)  mod  (2m-l) 

Thus,  (p,q)  are  defined  above  except  for  the  following  cases: 

(a)  if  n  even,   set  p  =  n-1        in  PE(n  -  1  -  k) 

set  p  =  n  -  1  -  k  in  PE(n  -  1) 

(b)  if  n  odd,  note  that  p  =  q  in  one  PE.   This  fact  will  be 
taken  into  consideration  later  on  in  the  program. 

As  the  pattern  of  pairs  is  constant  for  each  sweep  this  calculation  is 
only  done  once  in  the  program  and  saved  for  subsequent  usage. 

At  label  SWEEP,  all  necessary  preparations  are  made  for  another  2m  -  1 
t  r ans  forma  t  ions . 

(2)   ANGLE  -  Compute  sines  and  cosines  of  the  transformation  matrix. 

Input  to  tnis  routine  is  the  matrix  A  and  the  pairs  (p,q)  deter- 
mined in  (1)  above. 

Element    a^        is  brought  to  PE  q.   This  is  accomplished  by 
a  right  route  of  (q-p)  for  p  <  q  or  by  a  left  route  of  (p-q)  for 
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The  sines  and  cosines  are  computed  using  the  formulas  in 
Appendix  A  and  stored  in  two  rows  of  PE  memory  (SIN  ,  COS).   If 

■t-Vi 

n   is   odd, the   q        PE  has   p   =   q   and    in   this   PE,    cosine  and    sine   are 
set   to   1.0  and   0.0   respectively. 


(3)      MULTPL   -   Compute  B   =  A       R    .      Let   n  =  4,    for   k  =   1    the 

(k) 
ordered    pairs    are    (0,1),    (2,3).      Let  R..  be   as    defined    in 

Section  3.2.1  with  modifications   to   that  definition  as  noted 

in   (1)    of    this   section.      We  wish   to   compute: 


B  -  ViRk 


(k-1)     (k-1)    (k-1)  (k-1) 

00  01  02  03 

(k-1)      (k-1)      (k-1)      (k-1) 


l01 


11 


l12 


13 


(k-1)      (k-1)      (k-1)      (k-1) 
02  a12  22  23 

(k-1)      (k-1)      (k-1)      (k-1) 


03 


13 


23 


33 


Roo 

Roi 

0 

0 

R01 

Rll 

0 

0 

0 

0 

R22 

R23 

0 

0 

"R23 

R33 

We   note   that   column   1   of    B   =    [Rnn   x.col   1   of   A]+[ (-R      )    x   col   2    of  A] 


00 
column   2   of   B   =    [RQ,    x  col   1   of   A]+[R 

etc . 


or 


x  col  2  of  A] 


Rather  than  working  with  columns  it  is  preferable  to  work  with  rows, 
and  without  loss  of  generality  the  transformation  matrix  R  above 
may  now  be  considered  as  Rc. 


Then  B 
.    .    B 


(AR)L   = 


RtAt   =  r'a 


(since  A   is    symmetric) 


'V 


A, 


k-1 


Roo  Roi  °    ° 

-Rw   R        0     0 


0     R22R23 
0   -R23   RH 


(k-1) 

(k-1) 

(k-1) 

(k-1) 

aoo 

aoi 

a02 

a03 

(k-1) 

(k-1) 

(k-1) 

(k-1) 

aoi 

all 

a12 

a13 

(k-1) 

(k-1) 

(k-1) 

(k-1) 

a02 

a12 

a22 

a23 

(k-1) 

(k-1) 

(k-1) 

(k-1) 

a03 

a13 

a23 

a33 
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In  PE  memory  we  have: 

row  COS: 

row  SIN: 
row  PROW(p); 
row  PEN(q): 

Row      A    (0) 

A  (1) 
A  (2) 
A    (3) 


PE  0      PE    1      PE    2      PE    3 


Roo 

Ell 

*22 

R33 

R01 

"R01 

R23 

"R23 

1 

0 

3 

2 

0 

1 

2 

h    1 

1 

i 

1 
i 
I: 

Aoo 

Aoi 

A02 

I 
1 

A03 

1: 

Aoi 

All 

A12 

A13 

: 

A02 

A12 

A22 

A23 

: 

A03 

A13 

A23 

A33 

The   computation  of   B      is    simply  done   as   follows: 

Row  q   of   B     =  element   q  of   row   "COS"  x  row  q  of  A 

+   element   q   of   row   "SIN"   x  row  p  of   A. 

As  each  row  of  B  is  computed  it  is  skewed  in  preparation  for 
realignment  to  yield  B.  The  main  routine  needs  only  to  shift 
row   i  left  by   i    (1=0 ,1 , .  . .  n-1)    to   achieve   the  desired  matrix: 


b00  b01  b02  b03n 

b13  b10  bll  b12 

b22  b23  b20  b21 

b31  b32  b33  b30 


b00  b01  b02  b03n 

b10  bll  b12  b13 

b20  b21  b22  b23 

Lb30  b31  b32  b33 


B  skewed 
(output   from  MULTPL) 


B  re-aligned 
(in  nain  routine) 


(4)      Compute  R        B   =   A   .      The  routine  MULTPL    is    employed    to 
determine   the   new  A   .      Naturally   the  skewing   logic    in  MULTPL, 
described    in    (3)    is   bypassed. 
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(5)  Compute  R,     E      .    =   E    .      The   eigenvector  matrix   is    updated 
using   routine  MULTPL   once  more,    and   of    course  bypassing    the  skew- 
ing  logic.       It    is   preferable   to    pre-multiply ,    as    a  result    the 

rows   and    columns   of    E,    will   contain   the  left   and  right    eigenvectors 
corresponding   to   the  diagonal   elements  of  A,  . 

(6)  Convergence      (see   Section  4.2). 

(7)  GERSH  -  Finally   the  bounds   of    the   eigenvalues   are  determined 
in  routine  GERSH  by  Gerschgorin  discs.      The   eigenvalue   is    the 
center   of    the  disc,    and    the  bound  on   eigenvalue   i   is    the   sum  of 
the   off-diagonal    elements   of   row   i   in   the  diagonal  matrix  A,  . 
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h . 5     Program  Utilization 

The  usage  of  JACOBI   is   enabled  by  the   call   statement  below: 

CALL  JACOBI    (<  the    {,.    A  }  matrix  >, 

L  diagonal-1 

<  temporary  matrix  1   >, 

<  temporary  matrix  2   >, 

<  eigenvector  matrix  >, 

<  bounds  on  eigenvalues   >, 

<  order  of  matrix  >, 

<  sweep   count   > ) ; 

All  parameters    are  passed  as   addresses  whose   contents    are   described  as 
follows : 

1.  <  the    /  n .   A       .,  \   matrix  >   -  As   input   this   is  the  original  symmetric 

'•diagonal-' 

matrix  to  be   diagonalized.      If  the   user   desires  to   display  the   original 
matrix  or  retain   its    contents,   he   should  make  necessary  preparations  before 
the    call  to  JACOBI.      The   diagonal  matrix  replaces  the  A  matrix  and  is   output 
via  this   calling  parameter. 

2.  <  temporary  matrix  1   >  -   a  temporary  matrix  to  enable  matrix  multi- 
plication  for  JACOBI  which   is   available  to  the   user  after  exiting  JACOBI. 

3.  <  temporary  matrix  2   >   -   a  temporary  matrix  to   contain  the   pattern 
describing  the   elements   to  be   annihilated  in   a  given   sweep.      This  matrix  is 
also  available   for  usage   upon  exiting  the   routine. 

k.    <   eigenvector  matrix  >   -  the  matrix  of  eigenvectors   computed  by 
JACOBI  whose  rows    contain  the  right   eigenvectors   and  columns   the  left   eigen- 
vectors . 

5.  <  bounds   on   eigenvalues    >   -   a  vector  whose  value   in  PE  i   gives  the 

bound  on   X . . 

i 

6.  <   order  of  matrix  >   -  the  order  of  the  matrix,    an   integer  <   6k. 

7.  <   sweep  count    >   -  the   number  of  sweeps   to   achieve   convergence,    an 
integer  value. 

In   accordance  with   ILLIAC   IV  subroutine  linkage  standards,   the   contents 
of  ACARS  0   and  1,    as  well   as   $D32-$D63  are   saved.      The  user  is    advised  not   to 
use  $D0-$D31   since  they  will  be  overwritten. 
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PARAMETER 


STORAGE       STORAGE  CONTENTS 

REQUIRED      MODE       RELOCATABLE    DESTROYED 


<  the  ■[,.  A   ,  }  matrix  >     N  rows        Straight        Yes 
1  diagonal J  to 


Yes 


<  temporary  matrix  1  > 


N  rows 


Yes 


Yes 


<  temporary  matrix  2  >        N  rows 


Yes 


Yes 


<  eigenvector  matrix  >        N  rows        Straight 


Yes 


Yes 


<  bounds   on  eigenvalues   >  1   row 


Straight  Yes 


Yes 


<  order  of  matrix  > 


1  word 


Yes 


No 


<  sweep  count  > 


1  word 


Yes 


Yes 


where  N  =  <  order  of  matrix  > 
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5.0     Test   Results 


5.1     System  of  Even  Order 


See   [2]    p.    55.      The  matrix   to   be  diagonalized    is 


1  1 

1  1 
4  2 

2  4 


A2   = 

A3   = 
A,    = 


10.0 
5.0 
2.0 
1.0 


—     -> 
2 

-1 

""  0 

xl = 

2 

X2   = 

-1 

X3    = 

0 

1 

2 

-1 

1 

L.      J 

L-2- 

• 

1 

L 

x,     = 


Output  from  JACOBI  (note  the  eigenvector  matrix  could  be 
scaled  to  yield  the  above  result) '. 


*     Due  to  the  slow  speed  of  the  ILLIAC  IV  Simulator  (about  10   times  slower 
than  the  ILLIAC  IV),  we  tested  only  matrices  of  small  order. 
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5.2  System  of  Odd  Order 


See  [2]  pp.  58-59  the  matrix  to  be  diagonalized  is 


5 

4 

3 

2 

1 

4 

6 

0 

4 

3 

3 

0 

7 

6 

5 

2 

4 

6 

8 

7 

1 

3 

5 

7 

9 

A   =  22.40687532 
A2  =  7.513724155 
A  =  4.848950120 
A4  =  1.327045605 
A,  =  1.096595181 
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5.3     Comparison  with  Existing  Jacobi   Algorithm 

W.    Bernhard's   ILLIAC   IV   routine  EIGEN   [1]    is   essentially  the 
algorithm  discussed  briefly   in    Section  3.1.      A  comparison  run  was   per- 
formed  to    insure  the  two   algorithms   produced   compatible  results.      See 
[1]    pp.   127-129. 

Output  from  JACOBI: 
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APPENDIX  A  Determination  of  the  Rotation  Angle 


(k) 

The  orthogonal  matrix  R(p,q,a        ),   differs   from   the   identity 

matrix  by  a   2  x  2  diagonal   submatrix  whose  elements  are 

(A.l)  R        =   R        =   cos    a  R        =   -R        =   sin   a 

pp       qq  pq  pq         qp  pq 

(k) 
where  p     <   q.      In  order   to   eliminate  the  off-diagonal   element  a     "   , 

pq 

the  angle  a   is  chosen  such  that 

pq 

2a(k) 

<A'2>  tan2apq)=  (k)    ^    (k) 

rn  a  -  a 

pp         qq 

(k)  (k) 

in  which  a  is  restricted  by     a  <     /4.      Let 

pq  '  pq      = 

t     =      |J.(»|    ,        x     -      |a(k>-a(k>|    ,     , .    -(t?  +  x?)   V2   ; 
k  pq  k  '    pp  qq    '  k  k  k 


then 

(A.  3)  2      (k), 

cos     a 

pq 


1    /  m     "A  •    2      (k)  1    (_  Xk  "\ 

—   [1  +     — I  ;      sin   •  a  =—  |1-     —    /• 

2  V      yk/  pq       2  v      yk  / 


i  (k)  i  (k) 

Since     a   '  i      <      /4,  then  cos   a       will   always  be   taken  positive  and 

pq1    =  pq 

(k)       .„    .         c  .  .  ro    (k)  ,,    (k)         (kK, 

sm  a  will  be  of  the  same  sign  as    I  2a        /(a        -  a        )  J . 

pq  pq       pp       qq 
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APPENDIX  B 


Examples  of  Elimination  Scheme  [3] 


Let  n  =  8  and  k  =  2,  then  the  pairs  (p,q)  are  given  by  {(2,3);  (1,4) 
(7,5);  (8,6)}  and  R  of  the  form 


,(2). 


,<2) 
14 


R(2)    R(2) 

_.p(2)    r(2) 
K23      R33 


-R(2). 
R14 


(2) 
'44 


R(2) R(2) 

R55  ,57 


,(2)_ 


!2)_ 


Kl 


66 


G8 


"77 

(2) R 

68  R 


(2) 
33 


while  for   k  =   7   the  pairs    (p,  q)    are   {(8,1);    (7,2);    (6,3);    (5,4)} 
and  R     is   the  form 


H(7). 

!     ,u 


R(7) 


33 


,<7) 
.18 


R(7) 
K27 


?(7) 

*36 


R(7)  p(D 

R44  R45 

-R(7'  R(7) 

R45  *$S 


(7) 
'36 


(7) 
^66 


(7) 
l27   " 


R-. 


(7) 


'10 


(7) 
l88 
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If  the  order  of  the  matrix  is  odd,  say  n  =  7,  then  for  k  =  3  the 
pairs  (p,q)  are  given  by  {(1,2);  (7,3);  (6,4)}  and  R_  is  of  the  form 


R(3)  R(3) 

Kll  R12 

.R(3)  R(3) 

™12  '<22 


,(3). 
,33 


,  (3) 
,44 


^(3) 
^46 


R37~ 


,(3) 
137 


,(3) 
.46 


,*(3) 


!(3) 
'77 
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APPENDIX  C  A  Typical  Calling  Program 
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?USER 
7C0MP 
?BCL 


DEFINE 


I  *** 

MATAt 
MATB' 
MATP» 
E.I6V» 

&ERSH 
N« 

S*P» 

MtSSO 

MESS2 
MESS3 

^Essa 
mEss5 

START 
% 


•cacc 

ILE  H 

Card 
r£g 
fil 
pri 

CLC 

CAD 

CAD 
CSH 
CAD 
CCB 

CLC 
CAD 
CSH 
CCB 
DIS 

CAD 

CRH 
CAD 
CRO 
T*E 

BLK 
BL* 
BLK 
BL* 
'  BL* 
WDS 
WOS 
DAT 
DAT 
DAT 

DAT 
DAT 
DAT 
WDS 
FIL 


IGEN 

CD/A5K0/jArnBnRIVER  WITH  ASK  LIBRARY 


IN 

L 

NTMT*' 

Den 

Ld) 

D(l  ) 

CP 

(0)1 

D<0> 

L(0) 

(0) 

PLAYR 

<2) 

D(I) 

TR(i> 

D<1  ) 

TL<1> 
PM(o) 


12*| 


*Csl 
8D41I 


24 


*C3I 
151 

*D4l  i 

24) 

15, 

*Cl»"l6| 

■  6a» 

SC2| 

24, 


*C2> 
24  J 
•-gift1 

###f#####ENn  D£FINE#«##*#######*###*### 
INE  MMAxsl6f#) 


MMAXl 
MMftXl 

NMa^I 

NMaX  » 

II 
1» 
1> 

OtO#0» 
"#f#THE  ORIGINAL 
Mf##TME  DlAG^NAL 


MA Tp  I  * ###"> 
MATRlX#f #"| 
«#*#THE  EIGENVECTOR  MATRIX") 
w###MUMBfp  Or  SWEEPS  REQUIRED***"' 
"###  Pe  I  CONTAINS  BOUND  ON  F\<AL(I) 
0| 
L? 


#»#  ") 


SETE  E,TR,-F»  sET^l  E.AND.E' 

LIT(O)    l.N.N)   *  READ  DIMENSIHN  pF 

I  NpiiT  $c0.1j       *    SySTFm    FROM    INpnT 


LlT(O) 

STL(O) 

SL1T(0) 

LDAD(0) 

CSUB(O) 

STL(O) 

LIT(O) 

CRDTR(O) 

C*DD(0) 

CROTL(O) 


el  I 

SD40I 
■  N| 

SCO) 
SDftOJ 
SDttl  I 


*  FlXfD  PT  1  FOR  PRNTMTV  OEF 


%    N-l  FOR  PRNTMTV  DEF 


1. mata»mata* 

24| 

fD4l >    *  SET  UP  AcAR  0  FOR 

24)      *  USED  TO  READ  MATA 


INPUT  INSTRN 


00000100 

00000200 

00000300 
OOOOOftOO 

00000500 
00000600 

00000700 
00000800 
00000900 

ooooiooo 
ooooiioo 

00001200 
00001300 
00001400 
00001500 
00001600 

00001700 
00001800 
00001900 
00002000 

00002100 
00002200 
00002300 
00002400  I 

00002500  !:: 

00002600  I! 
00002700  j 
000028001 
00002900  I: 
00003000  Ij 

oooo3ioo  n 

00003200  ,j 
00003300 
00003400  I 
00003^00 

00003600  *j 

00003700  |j 

00003800 

00003900 

00004000 

00004100 

00004?00 

00004300 

00004400 

00004500 
00004600 
00004700 
00004800 
00004900 
00005000 
00005100 
00005200 
00005300 
00005400 
00005500 

00005600 
00005700 
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LlT(l) 

CAODM 


CPOTL(l) 

% 
RDINPTICLC(2)I 

CADD^?5 
CSHPC2) 
STA 

INPUT 
CROTR(O) 
LIT(22 
CADD(6^ 

CROTLCOI 

CADD<0) 

TXLTH(n 

LIT(0) 

DISPLAY 
CLC(3)I 
SLlT(3) 

*  ******* 


0»1 tOl 

%on\t 

?4  j 


t    ACAR  1  LOOP  INDEX 

*  FRO*  0  TO  N-l 

*  CLEAR  RSA  FOR  zFRn-FILLJNG  HTx 


SCO* 
61 

0(2)1 
*COM  I 
241 

■  641 

tC2| 
2a; 
to?* 
,RDInpTi 


ACAR2  USED  FOR  PE  ROW  A*0R 
CU  ADDR  FROM  ACARO 
PE  ADDR  TO  $C2 
ZERO-flLL  RDw  nF  mATPIx 

READ  MATRIX  FROM  INPUT 
RUMP  ACARO 

to  prepare 

r0R  NExT 

ROW  OF  MATRIX 

INPUT 


ImEss1"1  *MEssO| 
»C0.3?> 

«MATA,  pRlNTMTxi        *    ThE    ORlfilNAL     MTx 


%    THE    EIGENVALUES 


call  jacobi<mata»"atb.matp»figv.gfpsh»n.swp)j 

******* 

LIT  (0)  1  .'1ESS2-1  »MESSU 

DlSpLAY  &C0.32I 

SL1T(3)  -MATAI  P^INTMTXI 

tlT(O)  1 .MESS5-1 .MESSAI 

DISPLAY  tCo»3?» 

LlT(O)  1,gEpSH»p,ERSH! 

CROTR(0^  2ai 

cadd(O)  soan 

CRUTLCCO  24» 

oISpLAyR  8C0»16J  *  BOllNpS  ON  EIgENvALS 


?END 


LIT(O)  1»mES*3-1 .MESS2I 

DISPLAY  *C0»3?> 

SLIT(3)  -ElGVj  P^INTMTXj 

LlT(O)  1,mEss"-1 .mESS3| 

DISPLAY  *C0»3?» 

LlT(O)  I.mESSO-1 .$wp« 

DISPLAY  SCO • 1 6* 

HALTl 

FNC  START. 


*  THE  FlGENVECTnRS 


*  ITERATION  COU^T 


00005fl< 
000059C 
000060( 
00006K 

00006?( 
000063C 
000064C 

000065C 
000066C 

000067< 

000068« 

000069* 
000070C 

000071C 

onoo7?( 

000073C 
000074< 
0000751 

0000761 

000077(1 

0000781 

0000791 

000080! 

000081 

00008?| 

000083< 

00008/11 

OOOO85J 

000086* 

0000871 

000088 

000089 

0000901 

000091 

000092 

000093 

0000941 

000095 

000096 

000097 

000098 

000099 

000100 
000101 
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APPENDIX  D  JACOBI  Listing  (ASK) 
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?USER«C ACt  IGEN 

7C0MPILF  MCD/ASK0/NFW.|AC1RI  WITH  ASK 
?DATA 


LIBRARY 


.ZERO1 
,ONE  I 
»TMQ  i 


rEgIN 

FILL 

pQu 

eQu 

EQU 


•SIXT^IFQU 
.N  »  pQU 
•Nl  I  EQU 
•H  I  ECU 
•M2  I  pQU 
.M?M1»  pQU 
»MaM2»  fQU 
,SvLDOp|E«U 

****counteps» 

,AmT«   FqU 

.bpunOifqu 

.EvNOOD|EOU 

tlCNTl  EQU 

.KC0Nv,EQl' 
.RATlOiFQU 


1  2  ^  > 

$00  » 
$Oll 

sdzi 

JD3| 
$051 
8  06  J 

$071 
$D^  I 
SD9» 
SDiO 

$otl 

TNDIrAT 
SD1  2 
*Dl3 
$0l4 
*  Dl  5 
SOl  6 
$0)7 

JDtP 


dimension 
n  minus  i 

2*M 
2*M-l 


OF  SYSTEM 


tRnUTJ  FQU 

****S^ITCHES 

,ADS*T:eQU  $0i9 

•  RswUhiEou  $q?o 

i***SAVE  PATTERNS 

,PEMOL)EiE«U  SD21 

•SaVI I  FQU  J022 

•  SAy/EliEQU  $D?3 

•svmouifqu  $d?<» 

t***SAVE  REGISTERS  0 

. A0R1 •  Fyu  $o?5 

,Adr2|  EqU  fr)26 

•INDE*IEQU  *V?7 

•INL00P»EQU  $0?B 

.RETURiEQU  S029 

.SAyEOjEQU  Jo30 

tSAVEltEQU  *0  3 1 

.SAV3»  E^U  1032 

.SaVE3ieOU  JO  3  3 

.AnRA|  EQU  to!* 

.AORBI  FQU  S|)  35 

.ApRP*  EQU  $036 

.AOREl  EQU  $037 

•ApRGI  EQU  fr)3fl 

.ADRII  EQU  $039 

*  ****  NOTE- 
A 

COS  A l   PL*  II 

SlNA •   PL K  II 

pRnwj   ri*  1  > 

TEMPI »  RL*  1  I 

TEMP2I  RLK  ii 

TEMPPI  RLK  11 


0  TO  N-l 


nRS 


%    LOOp  pROM 

AND  LIMITS 

*  ROUTING  DISTANCE 

INDICATOR  TO  EN  A  RLE  C0V VERGpNCF 
EvEN/ppD  TNoIcATOR  pOR  N 

iteration  counter 

convergence  factor 

ratio  qf  sum  of  sops  or  qff.dtag 

Squares  Op  nlAr,  eifmFnts 

64"N 


TFST 


TO  St«M  OF 


*  aO  IF  OFF-OIAC, 
%    ROUTING  SWITCH 


SUM,  S1  IP  OTAr'  SUM 
(0«RT.1«LFT> 


*  ENABLING  PATTFRN  FOR  Pf-S  0  TO  N-1 


ADDRESSES 


AOOP  OF  MtX  TO  rf  mui.tplipp 
AooR  Op  MJx  cnNTAINlNG  RESni  T 
INDEX  USED  IN  (P.O)  SETUP 


OF  MULT 


return  addr  to 

SAvE  REGISTERS 


CALL  ING  PROG 


Address 
aodr  of 
addr 

ADDR 

AooR 

ADDR 


OF 
OF 
OF 
OF 


PROG  CURRENTLY  FXpEcTS 
1^  ORnEP 


Op  ThE  INpiiT  mtX 

MTX  CONTAINING  TEMP  RF SUITS 

MTX  CONTAINING  PATP*  (P»0) 

EIGENVECTOR  Mjv 

ROW  TO  CONTAIN  ROi)NnS  ON  FyAL 

sweep  count  fo»  output 
cosa,sina,  *  prow  to  re  storfd 


%  COSINFS  OF  TRANSFORM 

*  SINES    OF  TRANSFORM 

x  row  index  of  element 

*  TFMP  STORAGE 


MTX 
MTX 

T^  BF 


ANNTHlLATEn 


00000 
00000; 
00000* 
00000* 

00000* 
00000* 
000007 
OOOOOf 
00000 

ooooic 

00001 
000012 

00001" 

00001* 

ooooi 

00001* 

000017 
00001' 
OOOOlQ 
000020 

00002 
00002 
000023 
00002 

00002 

000026 
000027 
00002 

00002 

000030 

00003 

00003 

0O0033 

00003* 

00003s 

00003* 

000037 
00003' 
00003< 
00004C 

OOOOM 
0000*2 
0000*3 
00004A 

00004 

OOOOAA 

000047 

00004* 

00004< 

000050 

00005 

00005 

00005 

00005a; 

000055 

00005 

00005 


PtNl 


mil 


I  '1 


j 

* 
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en  t 


DATA  0»1»2.3»4.5.*»7»8»9»10»11»12.13.14.15»16.17»18»19.20, 

2 1.  22.  23.  24 »?5.  26 .27, 28 .29. 3 0.3 1.32*3  3. 34, 35 .36. 37. 38, 39 

41. 42»a3.44»a5»46»47, 48.49. 50.5l»52»5*»54»55.s6. 57.5a, 59 


ummyi  wos 
dbsaviwds 


01 


60.61 .62.63* 
MLL  l 6  I 


101 


*  SAVE  AREA  FOR  SD32-SD30  ANn  $C0-«C1 


NGLE'  FILL' 


ROUTINE.  ANGLE  COMPUTES  THE  SINES  AND  COSINES  0^  THr  TRANSFRORMATI 

THE.  COSINES  AND  SINES  ARE  STORED  IN  TWr  ROWS  Or  ^    MEM--COSA  »  SI 

and  are  deterged  8y  the  values  pen  (0>  and  p*ow  (p5  in  each  pe 

»I»E,  A(pVQ)»  A(4,n).  AND  A(Q,Q)  ARE  OFTERmINF*  ry  P»0. 


STLO) 

LD* 

LDU(O) 
LDR 
IDA 

SBM 

CLC(2)I 

IAG 

SETC(l) 

STL(l) 

SAP? 

IME 

StTC(3) 
STL(3) 

SKIP 


•SAVE3*   «  RFTURM  ADDR 

PEnI 

.AnRA» 

*0(0)  » 

PEN  j 

PRQWj 


SC2« 
II 

.SAvFTI 

*C2* 

II 

.  sAyI  I 

.  B  i j  M  p  y  2  J 


%  MAIN  DIAGONAL  TO  RPR 

%  SET  I  BITS  WHERE  P<0 

%  RQA  a  ABS(Q-P) 

*  IF  N  IS  EVEN  .SAVI.O 

%  iF  n  is  odd  .savI  note*  pF  WhP*e  Psq 


BRING  ELEMENT  A(P,P)  To  PE  0 


SETC<3) 
ZE«T(3) 

LDEE1 
LDL(l) 

CAND(l) 
L0EE1 

RTL 

STM 

LDL(i) 

COMpC(i,j 

CAND(l) 
LDEE1 
L0LC3) 
CSUB(3) 

CSUB(3) 

RTL 

STK 

LD8 

L0L(3) 
LDEE1 
UHPX2iaLIT(2) 


*C2» 
II 

,8'iMpy2| 

*C3I 
.SaVFT* 
$C3  » 

sell 

0(2)1 

TEMPPI 
.SaVeH 

*C3I 

*Cl  » 
.Si*T4i 

$C2l 

»C2I 
0(3)1 

TEMPPi 
♦0(0) l 
.PEMOnEl 
$C3l 
«ll 


*  SEE  IF  RGAa$C2 

*  pe-s  on  where  rga*$c2 

%  turn  on  pe-s  where  rga«sc? 

x  p>q  do  right  r0utf  of  *c2 


AND  p<o 

TO  GET  A(P.P) 


*  Turn  on  pe-s  where  rga*sc2  and  p>q 

%  bring  rgr  rack  into  position  and 

%  do  left  route  of  *c2  (pight  64«*r?) 

%  TO  GET  A(P.P) 

*  BRING  BACK  MAIN  DIAGONAL 


T.iRN 
LOOP 


FIRST    N-l     PE-S    ON 
FROM    1    THRU    N-l 


LESSTA(2)    S05.ANG1 I 


COMPUTE  SINES  A\D  rOSTNES 


00005800 
.40*00005900 

.    00006000 

00006100 
00006200 

00006300 
00006400 

00006500 
00006600 
ON  M00006700 
NA  00006800 
00006900 
00007000 

00007100 
00007200 
00007300 
00007400 
00007500 
00007600 
00007700 
00007800 
00007900 
00008000 
00008100 
00008200 

00008300 
00008400 
00008500 

00008600 
00008700 
00008800 
00008900 
00009000 
00009100 
00009200 

0000*300 
00009400 

00009500 
00009600 

00009700 
00009800 
00009900 
00010000 

00010100 
00010200 

00010300 
00010400 

00010500 
00010600 
00010700 
00010800 
00010900 

oooiiooo 
oooiuoo 

00011200 
00011300 

oooiuoo 
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lITCl) 

ID* 

LDA 

MLRN 
*TA 
LDX 
LDA 

sbrn 

5TA 

EOR 

CLC(1)I 
IAL 
SETC(3) 

lqa 

SETC(n 

ZEHT(l) 

LOL(O) 

CAND(1  ) 
CE*OR(3 
OR 

SETC(l) 
ZERTCl) 
L0K2) 
C0R(2) 

STL(2) 

SETE 

SETE1 
N0ZER0|STL(3) 

LDA 

MLRN 

LDS 

LDA 

MLRN 

AORM 

CALL  SO 

LDS 

LDA 

sap; 

OVRN 
LITCl) 

I  AG 

SETC(2) 
ZEHTC2) 
LDEE1 

LDA 

LDL(2) 
LDLFl 
SVTEMPlsTA 
ADRN 
I  I  T  (  0  ) 
MLRN 
CALL  SQ 
.STA 

LUA 


•2,01 
PRQWi 
♦0(0} I 
SC  1  I 
TEmPH 
PE^I 
*0(01 I 
TEMPPi 

TEMP?l 
TEMPI f 


SCI 
II 

TEm 

Jl 

»Nf) 
#SA 

*co 

SCI 

TEM 


I 


P2« 


Zf  rOi 
VETI 
I 
I 

Pi! 


%    TFMPl«2*A(p,Q) 


%    TFMp2«A(QtQ)-A(p,p) 


%    SET  I  BITS  WHERE  TEMPI  *  ?  DIFFER  TN  MGN 


II 

•  NO 
.SA 
SCI 

•  SA 
•I. 
E.A 

•  Sy 

TEM 
Mi 
SA  J 
TEm 
SAj 
SSI 
RT64( 
SA  i 
TEm 


ZERO! 

VII 

I 

VT! 

AND«E { * 

ND.E!  * 

mPhi  % 

p?l 

% 


Pi! 
)l 


CHECK  IF  ANY  TFMP?«0.0 

J  RITS  ON  WHERE  TEMP2«0,0 

SET  BUS  IN  $cn  WHFRE  »<Q 

P<Q  AND  TEMP2-0.0 

CHANGE  MODF  STATUS  nr  THESF  PF-S 

SO*E  TEMP?aO  SEE  T F  CORRESPONDING 

TFMpUO  ANJD  SET  I  pITS 

BRING  TO  ACAR  1 

IF  NO  TEMPI  K    7    BOTH  «  0  JUMP 

NOTE  THESE  PE-s  IN  PATTFRN  TO  SET 

COS  &  SIN  ,1,0  AND  0,0  RESPECT. 

WILL  RE  USFD  SHORTLY 

Turn  off  pe-s  wherf  tf«pi»tfmp2«o 

TO  AVpIDE  ZERO  DlVlDF 
SAVE  pATTERN 

TFMP2*TEMP? 


TEMP1*TEMP1 

SORT ( TEMPI **2+TEMP?**?> 


*S» 
■1,01 

SCI  I 

II 

»svtempi 

SC?  i      * 

SCI  I 

.pempdEi 

SC2I 

TEMPPi    * 
SCi  I 
■  0,5, 

SCO! 
9TM(  )  !     % 
COsAl 

sen 


INSURE  ABS(COS)  ft  A«S(MN)  LE*  ltO 


SAVE  TEMP?/(ABnVE  ROOT> 


COSINE 


00011' 
0001K 
00  0117 
00011P 
00011<5 
00012C 

000121 
000122 
000123 
r>00124 
000125 
00012* 
000127 
00012R 
00012Q 

000130 
000131 

000132 
000133 
000134 
000135 

000136 
000137 
00013^ 
000139 
OOOHO 
000141 

00014? 
000143 
000144, 
000145 

000146; 
00014  7i 
00014R 
00014QJ 

000150* 

oooi5i; 

00015? 

000153: 

000154! 

0001551 
000156( 
000157. 
00015fli 

0001591 
000160( 
00016K 
00016?( 
000163' 
000164( 
OOOI65C 
0001 66^ 
000167<1 
00016HC 
00016QC 
0  0  0  1  7  0  ( 
00017K 


tl 
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SBHN 
MLHN 
CALL  SO 
STA 

LDL(3) 

LDLE1 
CHSA  I 
STA 

cLcm  i 

10L(3) 

LOtri 

L.DA 

STA 

I  ft  A 

STA 
ftSET'  LDL(3) 
I  Dt.Fl 
LftL( 3) 
EXCHL( 3 
********* 


TEMPP? 

sco{ 
rT64(  ) ;  %   Sinf 

SIna, 

•SvMnn J 

*C3  > 

sl\|A  s 
.SAVIl 
*C1  5 

c  o  s  a  i    *  on  s  s  i ,  o 

?C?j 

S  I  \l  A  j      %     SlNsO.O 

.pfTMnnE? 
$  C}} 

************************************ 


***** 

I  ####«##«#*M(jlT IPlY  routine  **##«****######* 


MU 
MU 


* 

h 

h 
MULTP 


RE 


MA 

MA 


LTIPLY  Tr 

LTfPLY  TS 
ELEMENT 

♦  ELLMENT 
SUIT  IN  R 


TX 
TX 


IN  AftR 
IN  AftR 


L  J  FILL  ? 
STL<3) 

I  DUO) 
LDX 
fULTj   SLIT(1) 
CADftt 1) 
LOAOC1 ) 
LftL(3) 
fAOft(3) 

LQA 

ML*N 

LOS 

rADD( l) 

LOAO<l ) 
C  A I J  D  (  1  ) 

LOADC1 ) 

LDL(  1  ) 

CAUD( 1) 

LnA 

K'LHN 
AftKN 

*   IF  This  is 

LftL(  1  ) 


ANsF'iRMATI  im  MATRly  BY  MATpTv  Im  ApQl  --pFshLt 
DHNf  as  FlLLOwS-- 

:PF(g:  n>"  RflW  rns  *  «Ow  «PFn:  PF  mtv  IN  &DR1 

:  PE  n  s  n^  Rftw  MN  *  ROw  :PRftW:  nF  mtX  jn  Aor1 

T  W  |  P  F  N !  nr  MTX  IN  AftR? 

CALL  1     CALL  ?     0U.L3 
l    mata      matr      etcw 

■>  MAT6        MATA        MAjn 


tRA^F V 

.SVLftPP'     *     LftftP    FROM    ft     Tft    N - 1 

pEn  1 

■  C  0  5  A  I 

*CO? 

*c? s  *   ens   (P.Q) 

.An»l  : 

SCnj 

0(3);  *    LOAft    R[JW     (PEN)     OF    MATX     TN    ApRl 

*C2* 
*A  ; 
.SlXT/iJ 

*c?;   *  st.vj(p.o) 

.  S  I  X  T  U  t 

SC.31      *  "VALUF  OF  PROW  TN  PF  n 

•  A  n  R  1  1 

*C3; 

0(D!      *  LH  AD  ROW  (  P  R  ft  W  >  OF  MTRV  IN  AftRl 
*C?  ! 

pCJST    MnL  T  I  DL  Y    (CALL1'     "    3  K  E  W    M  A  T  P  T  X  " 
,Atri  : 


TN  MATRIX 


0O0172O0 
00017300 

000l7<j00 
0OO17S0O 
00017600 

00017700 
00017*00 
00017900 
0  0  018  0  0  0 
0  0  0 1  8 1  0  n 
00018?00 
00018300 
0001 8^00 
00018SOO 
00018600 

0  0018  7  0  0 
00018ROO 

00018900 
00019000 

000191 00 
0  0  0199  0  0 
0  0  019  3  0  0 
00019400 
IN0O019S00 

00019600 
00019700 
00019800 
00019900 
000?0000 

0  0  02  010  0 
00020200 
00020300 
00020^00 
00020SOO 
00020600 
00020700 
00020800 
00020900 

00021000 
00021 100 
00021200 
00021300 
00021400 

00021«SOO 
OO0216OO 

00021700 
00021800 

00021900 

00022000 

000221 00 

0002??Oni 

00022301 

00022400 

00022500 

00022600 

0002270  0 

0O022H00 


-39- 


MULTOl 


LDL(2) 

eQlxf(i  ) 

CLf(3)  i 

FGLXT(O) 

STL(O)    .AmT| 

SLIT(3)   »RnijTr* 

EXCHL<3^  SIpR  * 
STA 

CLRA, 

LDA 


. AoR?| 
»D3^#MULT1 I 

*C3tMnLT^iX  NO  NEED  TO  SkFW  FIRST  ROw 

*  skew  matx  in  aoR2  to  pfparf  for  transpose 


mULTI  « 
CKACOJ 


****** 

*  #### 

%   $A 

ROuTEl 


LDL<3) 
STL(3) 
cLC(3n 

SLIT(3) 

EXCHLC3) 
SWAPj 

LDA 

SKIP 

CADDf ?) 
STA 

TXLTM(O) 

LDL(3) 

EXCHL<  3) 


*0(2)j 
*X, 


*  STORE  RGA  SKEWfD  IN  MATY  fig  A0R2 


sot ; 

.AMTj 
■  ROUTE"! 

fiB.i 

•CKACOI 
*C0* 
0(2)1  * 

•MuLTi 

•SAVE3! 
MfRI 


%    ROUTE  PE  INOICTES  1  RI^HT 
*  No  DIRECT  PATH  FOR  LD*  *A 

STORE  RGA  IN  MTX  |N  ADR? 


******************************************** 

###*######t####pOUTE  ROUTINE  FOR  N-64 ###### *##############* 
HAS  ELEMENTS  TO  BE  ROUTED,  .AMT  CONTAINS  ROUTING  DISTANCE 


fill; 

STL(O) 
STL(1  ) 
STL(3) 
LDL(O) 

PTL 
LDA 
LOS 


*   SEE 


.SAVFO! 

.5AVF1  ! 

.  S  A  V  3  l 
•  A  M  T I 

SA.Of ni! 

PEN} 

0^ 


IE  LEFT(*D3i«l) 
CLC(l); 

EQLXT(  1  1 
LDL(3) 


RIGHT  ($031*0)  ROUTF 


RlGHTt 
SiTIT* 


CADD(O) 

rsup(o) 

ISG 

SKIP 
LDL(3) 
ISL 
SETF 
SLTF  1 

tlKa; 

RTL 
LDA 

LDL(3) 

L0U1 

LDL(O) 

LOL(  1  ) 

L0L(3) 

FXCHL(3) 


SD20.RIGMT! 
.Nj 
.Nl> 
.SfXTA ! 

*co» 

»SETITJ 
.RQIITi     %    64-N 

*C0! 

I .  ANn.El 

E. AND. El 

0<3)! 
f)R  | 

.pFMonEi 

JC3I 

.saveoj 

.SAVF1  ! 

.SAV3I 

IIrR| 


0O022Q 
000230 

000231 
00023? 

000233 
00023A 

000235 

000236 
000237 
00023A 
000239 

00024O 
000241 

00024? 
000243 
00024A 
000245 
000246 

000247 
000248 
000249 
000250 
000251 
00025? 
000253 
000254 
00025S 
000256 
000257 
00025A 
000259 
000260 
000261 

000262 

000263 

00026a 

00026S 

00026A 

000267, 

000268 

000269 

000270 

000271 

000272 

000273 

00027A 

000275 

00027ft' 

000277i 

000278 

0O027Qi 

000280) 

000281' 

000?8?i 

0002831 

0002841 

00028SI 


i( 


II 


-1+0- 


* 

31*****************  +  **  +  ***************  +  * 

audit:  fill i 

STU3)    ,SaVf1I 


)ITI  FILL  I 

ADOIT  CALCULATES  THE  *i)MM  OF  THF  PFF"PI  AGON ALS  SQ 
DIVIDES  IT  qY  THE  SUM  nF  THF  DIAGONALS  SQUARED 


IIAPF  ANP 


C  U  C  (  1 )  | 

STL(l) 
LDL(0^ 
SETf 

SETE1 

rLRA, 

LPS 

LDUO) 

LDLH 
LDL(2) 

rAUI)(2) 

LP  A 

ML^N 
AOHN 
LOS 

TXEfM(O) 


»  A  r)  s  w  T  i 

•  Sv/LnnP; 
E  •  n  p  .  -  E  j 
E  .  a  N  r  ,  E  ; 

*A  t 

,PFMnnE» 

.A  PR A: 

«COJ 

n(  ?)t 
*A  t 

*A; 

»  A  f)  I  j 


%   End  rpw-siim 


SETF  E.lR.-r» 


SETE1 
LPL(O) 

LUS 
RTL 
t  PS 
AnRN 
CAQP(O) 

LDK1) 

E®LXF(oi 

LDL(l) 

FQLXT(  1  ) 

S1  A 

PLRA? 

LOS 

LDL(l) 

L  0 1  E  1 

LDX 

L0L<2> 

LOA 

|^LKM 

LDL(l) 

$TL<  15 

SKIP 

LPS 

LDA 

SB-RN 

cLCcon 

IHL 

SETrtO1* 
PNEST(O) 
^VK^ 

LDC(O) 
STL(O) 


E. A^P.E  j 
,'OnFj 

%  A  .         * 

*S,0(O) *  * 

*R  j 

SSj 

*  c  ft  I 

.  S  I  X  T  4  * 

ICl  ,  Anil  >*  END  LtWi-SUM 

.  A  0  S  k  T  ? 

*  P 1 ,  A  n  1 2  : 

TFMpp, 


SUM  lip  ThF  ElFmEmTs  IN  EAcm  p? 
USFP  FOR  PPTH  nFF-PTAGnNAl.  «;  AND 


DIAGONALS 


*  a: 

.PFMnnEj 

*  r.  1  ? 

PFviJ 
.  A  o  R  A  l 
*0(2)! 
%  A  ; 

.PnEi 
. AoSWT? 
»  A  r)  I  0  ! 

*  A; 

T  E  M  P  P  ! 
«.S| 

*Co  s 
IJ 

•  A  n  I  3  I 
*SS 

*A» 
.RATTH J 


*  pick  up  the  Ari.ii-.s 


X  SnRTRAcT  TmE  Af  I.  I  1 


%    CnNVERdENCF  FArTnR 


ooo; 
ooo! 

poo? 

ooo? 

0002 

0002 
000? 
0002 
0002 

000? 

ooo? 

000? 
000? 

0002 
0003 
0  00  3 
0003 
0003 
0003 
0003 
0003 
0003 
0003 

0003 
0003 
0003 
0003 

0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 

0003 
0003 
0003 
0003 

0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 
0003 


8600 
8700 

8R00 

8900 
9000 

9100 
9?00 
9300 
9400 

9soo 
9600 
97  00 
9*00 

990r 
0000 
01  00 
0?00 
0  3  0  0 
04  00 
050  0 
0  6  0  0 

0  7  0  0 
0800 

0900 
10  0  0 
1100- 
1?00 

1  300 

1  a  o  c 

1500 
1600 

1700 
1800 

1900 
?000 

2100 
2?00 
2300 
2400 

25  00 

2600 
2700 
?800 
2900 

3000 
31  00 
3200 
3300 
3400 
35  0  0 

3600 

3700 
3800 
3900 
4000 
4100 
4?00 


-1*1- 


ADI3I 

******* 
GERSH* 


******* 

*  THE 

%       OF    J 

*  LOCA 

*  Of  E 

j****** 

******* 

*  find 

*  OF  T 

******* 

&ERSS 


******* 
******* 

JACOfiH 


SA  | 


LDL(3) 

E*CHL<3) 

*** 

FILL  * 


STL(3) 
SETE 
SETE1 
CL«A, 

LOS 
*** 

rauius  n 

HE  SUM  n 
TEL)  IN  9 

ach  Rnw 
*** 

LDLf  O 

*** 

I  N*»  THF 

HE  COLIM 
*** 

LUL(O) 
C  C  b  (  0  ) 
LDEE1 
LDL(2) 

CAUDC2) 

IDA 
sApj 

ADRN 
IDS 

LDt(O) 

LDEE1 

LDL(O) 

STA 

LDL(3) 

FXCHL(3) 

*** 

******** 

ENTRY):  j 

CHwS 

SETE 

SETE1 

STL(3) 

r  L  c,  (  3 )  i 

SLlT(3, 

STURE(3) 

A  L  1  T  (  3  ) 

RTOREf 3) 

CLC( 3)  ) 

SLIT(3> 
I IT(O) 
STOREOi 
ALlT(3) 
TXfc FM(o> 
LITCOJ 
CLC( 3)| 


. SA VE  3> 
MfiRI 


.SAVE3I 
E.DP.-EI 

E.anh.Ej 


x    GERSH  FINDS  A  RQUNO  ON  THF  F  JGENVaUIEs 
X    I)SING  THE  gERShgORINg  "IS* 


*A» 


F  THE  DlS*  ACC.  TO  GERSHGORIN-S  TH"  consists 

F  THE  ABS.  VALUE  Or"  THF  OFF-OIAf,  ELEMENTS 

OW     J     FOR     ^HICH    AiI»I)«FlGENVAL#  I  .Ff     THF    ROWSIIN 


5  pOuND.   T 

.svlooP' 


CH  A(I»I)wFIGEnVAL.I. 
HF  pENTER  nr  TMF  DISK 


K  TS  TgF  FTGFNVAL,|F 


SUM  OF  THC  ROW  IS  EQUIVALENT  TO  FADING  THE  SUM 
S.  SINCE  THE  MATRIX  IS  SYMMETRIC 

.PfMnnEj 
*co  j 

.  Ar)RA  I 
SC  1  X 
0(2)1 


*s» 

SA  I 
»GeRS> 

,pemopE j 
*Coj 
.  A  n  R  r,  J 

0(0)  J 

.Savet; 
McR? 


?!  STORE  BOUMOS  ON  EIGENVALUES 


*********************************** 
FILL?     *  SAVF  JC0.RC1*  AD«  S032-S039 
64, 

E.or.-Ei 
e.and.ej 
•retupj 

adrsa v  +  fl» 
*Col 

1 1 

%C.  1  I 

aorsav' 

1»  7»0I 
%  0  3  2  (  0 )  I 

1 I 

»sa» 
1 • 4»n  i 


SM 


000343 

00034a 

000345 
000346 
000347 
000348 
000349 
000350 
000351 

000352 
000353 
00035a 

000355 

^00356 
0003571 

00035* 

000359' 

000360 

0003611 

00036?i 

000363 

000364i 

000365 

000366 

O00367i 

000368* 

000369! 

000370 

000371 

00037?i 

000373< 

000374* 

000375 

000376! 

000377) 

O0037«i 

000379( 

000380! 

000381 

0003B? 

000383 

000384 

0003851 

0003861 

0003871 

000388< 

0003891 

0003901 

000391 ( 

00039?! 

0003931 

000394( 

0003951 

000396( 

000397( 

00039fl( 

000399( 


-1+2- 


L0AD(2) 
CShH(3) 

STLC3) 

AL1K2) 

T*EFM(0> 

LOADC?) 

LOADC3) 

STL(3) 

A  L  1  T  (  2  ) 
STU(3) 


SC3J 

*i 

*  0  3  4  ( 0 ) ; 

»sai> 

*C3  t 
SC  3} 

VM  s 
1  > 

*  c  3J 

,  a  n  R 1 1 


initialize  constants 


Sf2    SFT     IM    CALI      STMT    FTIM    MA  IN    PRpr, 
AnDP    [IF    DTAG'TFMP*     PAI"    AMn    FTf,V    MfXS 


-X    DlMENsInN    HF    SvSTEm 

t    AHOR       IJF     SWFFP    COUNT    -    TO    r>F    nWTPlJT 
A,sjn    hThEr    pApAmFtFrs 


CLC (Q)' 

STL(0i  .^FRn* 

STL(O)  .RSwTTH; 

STL(O)  .IrNT} 

STL(O)  ,Kcr^'VJ 

LIT(O)  «1  J 

stu  0)  .omf  ; 

cAUD(O)  *co; 

$TL(0)  .  Twn» 

I  I  T  ( 0 )  »  a ! 

5TKO)  .-RQUNn* 

I  I  T  (  0  )  s  6  4  J 

STL(O)  ,SlXT/iJ 

CSUHCO)  ,m, 

S  T  L  (  0  )  .  R  rj  1 1 T  : 

l.f)L(0)  ,  MJ 

cSUb(O)  mi  ; 

STL(O)  .Nil  *  N  MINUS  1 

1. 1 T  (  0 )  =  0  . 1  .  n  1 

CADD(O)  ,Ni  ! 

cRUTL(O)  24; 

STL(O)  .SvLnnPl 

C  A  U  0 ( 0  )  SOll 

CSHR(O)  1  j 

STL(O)  .Ml     *    M=r     (N*1  )/?     1 

CSHI(0>  1? 

STKO)  ,M?I     *    2  +  M 

CSUP(O)  ,N| 

STL(O)  .  E  y  N  fi  n  n  1 

LOL(O)  .M?J 

rsuRCo  .n^F? 

STL(O)  .m?Mij 

rSHL(O)  is 

STL(O)  ,M4M?:        %   a*M-? 

HN     FIRST  M»l     nr-S 


%    =0     lF-    N    Ey/FN,     =1     U     Pp 


LOA 

i_ol(0) 

IAL 

SETf.CO^ 
SU(0} 


P  E  m  I 

•  N  I 

*Co : 

T  I 
,PrMnnFl 


oooaoono 
n  0  0  a  0 1  0  0 

ooo40?oo 

0  n  0  4  0  3  0  0 
00040400 

00040500 
00040^00 
00040700 

00040«00 
0  0  04  09  0  0 
000  4100') 
00041 1 00 
0O041P0O 

00041 300 
00041400 
0  0  04  15  0  0 
0004160O 
00041 700 
00041R00 

0  0  0  4  1  9  0  0 
00042000 
0004210^ 

00042200 
000  42 3  00 
00042400 
0004250O 
00042600 
0  0  0  4  27  0  0 
00042800 
0004290O 
00043000 

00043100 
0O043?OO 
0004330  0 
0004340  0 
00043500 

OO043600 
0O043700 
00043flOO 
OO043Q0O 

OO044OOO 
0  0  0  4  4  1  n  0 
00044?0  0 
00044300 

0004440U 
0004450  0 
00044600 
00044700 

00044*00 
00044O00 
00045000 
00045100 
00045?00 
00045300 
00045400 
0004550  0 
00045600 


-k3- 


L 
INITI 

C 
L 


EylNlT | S 

C 
T 
L 

L 
L 

L 

S 


DE-El 
ALIZE  E 

LHA» 
0L(0) 

DL(l) 
TA 

ADO(l  ) 

XLTM(0^ 

D* 

|T(0) 
DA 

OL(t) 
TA 


SET  [)P    PAIRS 


DL(1} 
DEE1 


*C0f 
IgEmvEcTnr  MATRly 

•SVLOnP*  *  LOOP  FROM  o  TO  N-l 
. AnREl 

"(hi 
•  onf  t 

.EvlNlT* 
Pen* 
*1  .01 

*cn  ? 

.  AnREj 
♦0(1  )» 

(P.Q)  Fo<*  ANNIHILATION 

.PEMonE! 


scif     *  turn  on  pe-s  o  to  n-i 

SET  UP  LOOP  COuNTFR  t0  Go  FRO*  1  TO  ?"-! 


SETuri 


KNQTMI 


IT(O) 
ADD(O) 

ROTL(O) 
ADD<0> 
DL(1) 
TL(1) 

QLXFA(0 

AUD(l) 
SUB(l) 
TL(1) 

SUB(1  ) 

SUB<1) 

DA 

BM 

DL(2) 

AL 

ETF 

ETE1 

PM 

DL(2) 

DtEl 

DL(2) 

A0D(2) 

TA 

LC<3) J 

QLXFAM 

SHL(2) 

AU0(2) 

SUR(2) 

DL(  3) 

TURE(?) 

AUD(2) 


X  LOOP  FROM  1  TO  ?M-1 


■0,1 ,0* 
,M2M1 I 
24j 
*0l* 

.  iNOFyJ 


)  $n7,KNnTMj 

.M?J 
.OnFI 

•  Inofyi 

*coj 

*co> 

*C  i  ; 

PENi 

•M?Ml | 

*C?I 

-I  ,  AfvO.E  t 

E, ANn.E  I 

SC2;      *  IF  P   r.FQ   2*M-1SUR  2*M"1 

.pEMonEi  *  Turn  first  n-i  pf-s  on 

.AnRPl    %    qFT  ApoR  Or  PATR  MATRlv 
*C0 I      X  GFT  I-7H  ROW  AND  STORF  P»0 
0(2)  I 

%     IF  N  IS  000  SKjP  TO  POnnNF 
)  SnMtBlJMROI 


X  6*M-3   RESET  wmFN  RtM 


%     iNOfX  -  2*K 

X  IMOFX  -  2*k  -  0 


At 

•  Nl  I 

.Nl  | 
SC3I 

*CoJ 


x  back  to  cu  aohr 

*  ♦  ( n-i  ) 

X  -K 

X  STORF  N-l  IN  Pf(N-J-K) 

X  Cli  AOOR  PROW  ♦  (N-1  ) 


000457 
00045* 

00045<; 

00046C 
000461 
000465 

00046 

000464 

00046^ 

00046f 

00046? 

00046f 

O0046<! 

00047C 

000471 

00047; 

000471 

000474 

00047? 

000476 

000477 
00047ft 
000479 
00048C 
000481 
00048? 
000483 
000484 
000485 
000486 
000487 
000488 
000489 
000490 
0  0  0  4  9  t 
00049? 

000493 
000494 
00049«5 
00049fi 

0  0  04  97 
000498 
000499 
000500 
000501 
00050? 

000503 
000504 
000505 
000506 
000507 
00050ft 
000509 

0O0510 
000511 

000512 

000513 


Lfl 


!( 


-44- 


c 
s 

bUMpOl  I 


% 

SWEEPS  »L 

r 
s 
I 

L 


SUB(3) 

TCRE(2) 

DL(1) 

XLTM(O) 

E^  LOO* 

DL(0) 

AUO(O) 

TL(O) 

PEE1 


sen? 
•Setup  J 


N-t-K 
STflPt 


N-l-K  IN  PE(N-t) 


-  r|NF  swc-fp=?m-i 

ERF  M=[(^l )/2]  . 

.In  NT i 

,omE; 

.  I  r  N  T  ' 

,PEMHnE; 


I  TF RATIONS 

N  =  n  t  m  F  n  %  i  n  n 


or  cyStfm 


*  turn  on  pf-s  n  to  N"M 


BEGIN  ITFRATTOn  |  nGir 


.OOP* 


POST 

ACTUA 

THEN 

C 

S 
E 
I 

r 
i 

L 

C 

L 

S 

TRANS'  I 

S 

r 

s 

F 
A 

r 

T 


IT(O) 
AUD(O) 

ROTL(O) 

A  0  D  (  0  1 
TL(Oi 
DL(2) 
A  0 l>  (  2  ) 

OA 
TA 

L  C  (  3 )  J 

L1T13) 
XCHL(3) 

Dt(2) 
TL(2) 
0L(2) 
TL(2) 

MULTIPLY 
LLy  Prt- 
TAKE  TrA 

LC( 3)! 

LIT<3) 

*  C  H  L  (  3  1 

DL(O) 

A  L;  D  (  0  ) 
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